%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Assemble the matrix that maps an charge distribution on C to a
% outgoing potential
% 
% kh       Helmholtz parameter.
% C        contour
% ind      indices of the points to consider on C
% zz       positions where the potential is evaluated
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function res = get_B(kh,C,ind,zz)

M = size(zz,2);
n = size(ind,2);

% dist1(i,j) = x1(i) - z1(j)
dist1 = (C(1,ind)' * ones(1,M)) - ones(n,1)*zz(1,:);
% dist2(i,j) = x2(i) - z2(j)
dist2 = (C(4,ind)' * ones(1,M)) - ones(n,1)*zz(2,:);

Y_dg1 = (C(2,ind)' * ones(1,M));
Y_dg2 = (C(5,ind)' * ones(1,M));

nn1 = ( Y_dg2./sqrt(Y_dg1.*Y_dg1 + Y_dg2.*Y_dg2));
nn2 = (-Y_dg1./sqrt(Y_dg1.*Y_dg1 + Y_dg2.*Y_dg2));

dd  = sqrt(dist1.^2 + dist2.^2);

K   = funcK(kh,nn1,nn2,dist1,dist2,dd);

res = K.*sqrt(Y_dg1.^2 + Y_dg2.^2);


return
